Chapter 9 Cokriging

9.1 Librerías

library(sp)
library(gstat)
library(sf)
library(rgdal)
library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:MASS':
## 
##     select

## The following object is masked from 'package:vcdExtra':
## 
##     summarise

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(Matrix)

9.2 Descripción de los datos

Cokriging para las variables \(NO2\), \(O3\), y \(NOX\). La variable de principal riesgo es ozono (\(O3\)), así que se usan las otras dos como covariables espaciales. Día 2020/01/16 A las 17 horas.

datos <- read.csv("2_COK_G_stat/Air_polution_cdmx_2020_01_16_17h.csv")
datos <- datos[c("Estacion",
               "X",
               "Y",
               "NO2",
               "O3",
               "NOX")]

pander::pander((datos))

Estacion	X	Y	NO2	O3	NOX
AJU	482901	2117907	NA	50	NA
AJM	478188	2130946	5	51	8
ATI	473346	2164689	20	70	20
CAM	482180	2152665	23	83	24
CCA	481502	2136931	6	46	8
CUA	469366	2141275	9	56	10
CUT	479189	2180751	13	75	14
FAC	474444	2154232	31	71	35
HGM	484020	2146380	25	81	31
IZT	487647	2143367	21	61	23
LLA	495842	2164872	14	81	15
LPR	487650	2160000	NA	64	NA
MER	487445	2147815	23	79	26
MGH	478716	2145543	NA	66	NA
MON	510196	2151776	9	80	10
MPA	498809	2123036	1	50	NA
NEZ	497038	2144394	8	62	9
PED	478557	2136817	6	56	6
SAG	496819	2159801	16	80	16
SFE	472393	2140390	9	66	10
TAH	498890	2128098	3	47	3
TLA	478535	2159383	26	64	30
TLI	481421	2167509	18	86	20
UAX	489113	2134517	NA	54	NA
UIZ	492241	2140751	7	45	7
VIF	489875	2173664	19	80	22
XAL	491355	2159031	27	80	27

9.3 Matrices de coregionalización.

9.3.1 Matriz definida positiva para el modelo Esférico.

mat1 <- cbind(c(30, 30, 30),
              c(30, 50, 30),
              c(30, 30, 35))
#matriz definida positiva "cercana"
mat1 <- data.frame(as.matrix(nearPD(mat1)$mat))
names(mat1) <- c("NO2", "O3", "NOX")
row.names(mat1) <- c("NO2", "O3", "NOX")
pander::pander(mat1)

	NO2	O3	NOX
NO2	30	30	30
O3	30	50	30
NOX	30	30	35

9.3.2 Matriz definida positiva para el modelo efecto Hueco.

mat2 <- cbind(c(13.02, 24.5, 18.739),
              c(24.58, 46.4, 35.36),
              c(18.73, 35.36, 26.95))
mat2 <- data.frame(as.matrix(nearPD(mat2)$mat))
names(mat2) <- c("NO2", "O3", "NOX")
row.names(mat2) <- c("NO2", "O3", "NOX")
pander::pander(mat2)

	NO2	O3	NOX
NO2	13.02	24.54	18.73
O3	24.54	46.4	35.36
NOX	18.73	35.36	26.96

9.4 Definición de objeto en gstat

9.4.1 Semivariogramas univariados

vgmno2 <- vgm(psill = mat1[1, 1],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 1],
                         model = "Hol",
                         range = 2294))

vgmo3 <- vgm(psill = mat1[2, 2],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[2, 2],
                         model = "Hol",
                         range = 2294))

vgmnox <- vgm(psill = mat1[3, 3],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[3, 3],
                         model = "Hol",
                         range = 2294))

9.4.2 Semivarogramas cruzados (Bivariados)

vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 2],
                         model = "Hol",
                         range = 2294))

vgmno2_nox <- vgm(psill = mat1[1, 3],
            model = "Sph",
            range = 6096,
            add.to = vgm(psill = mat2[1, 3],
                         model = "Hol",
                         range = 2294))

vgmno3_nox <- vgm(psill = mat1[2, 3],
                  model = "Sph",
                  range = 6096,
                  add.to = vgm(psill = mat2[2, 3],
                               model = "Hol",
                               range = 2294))

9.4.3 gstat

remove_na <- function(frame, vari_) {

    # Remove na from sp object

    datos1 <- frame

    bool <- !is.na(datos1@data[vari_])
    datos1@data <- datos1@data[bool, ]
    datos1@coords <- datos1@coords[bool, ]

    return(datos1)

}

coordinates(datos) <- ~ X + Y

g_st <- gstat(NULL,
              id = "NO2",
              formula = NO2 ~ X + Y,
              model = vgmno2,
              data = remove_na(datos, "NO2"))

g_st <- gstat(g_st,
              id = "O3",
              formula = O3 ~ Y,
              model = vgmo3,
              data = remove_na(datos, "O3"))

g_st <- gstat(g_st,
              id = "NOX",
              formula = NOX ~ Y,
              model = vgmnox,
              data = remove_na(datos, "NOX"))
#Cruzados


g_st <- gstat(g_st,
              id = c("NO2", "O3"),
              model = vgmno2_o3)

g_st <- gstat(g_st,
              id = c("NO2", "NOX"),
              model = vgmno2_nox)

g_st <- gstat(g_st,
              id = c("O3", "NOX"),
              model = vgmno3_nox)


pander::pander(do.call(rbind, g_st$model)[, 1:3])

	model	psill	range
NO2.1	Hol	13.02	2294
NO2.2	Sph	30	6096
O3.1	Hol	46.4	2294
O3.2	Sph	50	6096
NOX.1	Hol	26.96	2294
NOX.2	Sph	35	6096
NO2.O3.1	Hol	24.54	2294
NO2.O3.2	Sph	30	6096
NO2.NOX.1	Hol	18.73	2294
NO2.NOX.2	Sph	30	6096
O3.NOX.1	Hol	35.36	2294
O3.NOX.2	Sph	30	6096

9.4.4 Estimación del semivariograma

plot(variogram(g_st),
     model = g_st$model,
     pl = T,
     xlab = "Distancias",
     ylab = "Semivarianza")

9.4.5 Mapas de predicción de O3 con las covariables espaciales NO2 y NOX

prediction_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 100000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)

    pred <- paste(variable, ".pred", sep = "")

    plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
            geom_tile() +
            scale_fill_viridis_c() +
            theme_void()

    return(plot)

}


variance_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 10000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)

    var <- paste(variable, ".var", sep = "")

    plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
            geom_tile() +
            scale_fill_viridis_c(option = "inferno",
                                 direction = -1) +
            theme_void()

    return(plot)

}

cv_plot <- function(g_object, variable, map_path) {

    map <- readOGR(map_path)
    new <- sp::spsample(map, n = 10000, type = "regular")
    coordinates(new) ~ x1 + x2
    colnames(new@coords) <- c("X", "Y")

    predic <- predict(g_object, newdata = new)

    prediction <- data.frame(predic)
    pred <- paste(variable, ".pred", sep = "")
    var <- paste(variable, ".var", sep = "")
    aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
    aux[aux > 1] <- 1
    prediction["cv"] <- aux

    plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
            geom_tile() +
            scale_fill_viridis_c(option = "magma",
                                 direction = -1) +
            theme_void()

    return(plot)
}


pl1 <- prediction_plot(g_st, "O3",
                       "2_COK_G_stat/SP/mpiosutm.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]

pl2 <- variance_plot(g_st, "O3",
                     "2_COK_G_stat/SP/mpiosutm.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]

pl3 <- cv_plot(g_st, "O3",
               "2_COK_G_stat/SP/mpiosutm.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]

ggplotly(pl1)

ggplotly(pl2)

ggplotly(pl3)